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Summary 


The  dispersion  of  a  plume  from  a  stationary  or  moving  gas  source  into  an  environment  is  representative 
of  accidental  or  deliberate  release  of  gases  from  land-  or  aerial-based  vehicles,  release  of  biochemicals, 
or  the  release  of  odors  from  biological  systems.  Dispersion  can  be  modeled  by  advection-diffusion  PDEs 
with  spatially  and  time  varying  ambient  mean  velocity  and  eddy  diffusivities.  While  mobile  sensor-centric 
schemes  for  gas  source  detection  can  adequately  perform  for  dispersion  processes  with  constant  wind  and 
diffusivities,  they  are  not  applicable  to  processes  with  spatiotemporally  varying  wind  and/or  diffusivities. 
The  model-based  detection  scheme  incorporates  the  moving  sensor  location  into  the  state-space  formulation 
of  the  dispersion  process  and  thus,  can  perform  under  general  dispersion  conditions.  Moreover,  the  increase 
in  computational  power  and  efficiency  in  data  acquisition  hardware,  enable  the  implementation  of  the  model- 
based  detection  scheme  onto  a  Sensing  Aerial  Vehicle  (SAV).  The  development  of  a  model-based  detection 
scheme  constitutes  the  main  areas  of  this  research  effort. 

The  goals  of  research  performed  under  FA9550-09- 1-0469  were  to: 

•  Develop  a  model-based  estimation  that  provides  the  sensors  position  employed  to  find  the  proximity 
of  the  moving  source.  The  sensor  spatial  relocation  is  accomplished  by  means  of  an  SAV  that  has 
navigation  and  communication  capabilities.  The  emphasis  of  the  effort  is  on  arriving  at  a  model- 
based,  optimal  sensor  repositioning  during  the  search. 

•  Develop  a  finite-volume  computational  method  on  unstructured  grids,  that  provides  in  realtime  the 
solution  (process-state  estimate)  of  the  2D  advection-diffusion  PDE  with  variable  diffusivities  and 
ambient  wind.  The  unstructured  grid  allows  modeling  of  complex  boundaries  (e.g.  objects)  and 
will  be  adapted  with  local  refinement  and  coarsening  during  the  process-state  estimation,  in  order  to 
improve  accuracy  and  efficiency. 

The  objectives  of  research  performed  under  FA95 50-09- 1-0469  were  to: 

•  Develop  a  new  theoretical  and  approximating  framework  for  the  simultaneous  estimation  of  the  source 
location  and  the  process  state  associated  with  a  stationary  or  moving  gas  source  that  emits  gases  in  an 
ambient  environment. 

•  Develop  a  performance-based  guidance  of  the  sensing  aerial  vehicle. 
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•  Develop  and  implement  a  computational  method  that  provides  in  real-time  the  solution  (process- 
state  estimate)  of  the  generalized  PDE  that  models  the  2D  advection-diffusion  process  with  variable 
diffusivities  and  ambient  wind. 

•  Theoretical  implementation  of  the  proposed  scheme  is  accomplished  via  an  SAV  that  provides  the 
spatial  relocation  of  the  sensor  and  communicates  with  a  land-based  station. 

The  Mathematical  and  Computational  Accomplishments  of  the  estimation  approach  are: 

1.  Examined  various  guidance  schemes  for  the  sensing  aerial  vehicle  carrying  the  sensors. 

2.  Examined  both  kinematic  and  dynamic  equations  of  motion  for  the  guidance  of  the  sensing  aerial 
vehicle. 

The  Theoretical  Accomplishments  of  the  estimation  approach  are: 

1.  Developed  a  model-based  estimation  scheme  that  provides  in  real  time  the  estimate  of  the  concentra¬ 
tion  due  to  the  gaseous  release. 

2.  Provided  an  estimate  of  the  proximity  of  the  moving  gaseous  source  via  the  motion  of  the  mobile 
sensing  aerial  vehicle. 

3.  Developed  a  vehicle  guidance  scheme  using  Lyapunov  redesign  methods. 

4.  Incorporated  the  motion  of  the  sensing  aerial  vehicle  into  the  estimation  equations. 

5.  Dictated  the  motion  of  the  sensing  aerial  vehicle  by  the  performance  of  the  estimation  scheme 

The  Computational  Accomplishments  of  the  computational  method  are: 

1.  Developed  a  multi  grid,  multi  step  finite  volume  method  with  upwind  and  flux  limiting. 

2.  Developed  grid  adaptation  with  local  refinement  and  coarsening  during  the  process  state  estimation. 
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Computational  Base 
Station  with  wireless 
communication 

Figure  1:  Overview  of  model-based  estimation  scheme  with  sensing  aerial  vehicles. 

1  Overview 

The  goal  of  the  work  is  to  provide,  in  real  time,  an  estimate  of  the  gas  concentration  associated  with  an 
emitting  moving  source.  Minimizing  the  damage  from  a  toxic  release  must  address  not  only  the  source 
location  or  its  proximity,  but  also  the  contaminated  material  that  has  already  been  released  (i.e.  estimate  the 
concentration  field).  The  strategy  applied  in  this  research  effort  is  not  to  reposition  the  sensors  to  areas  of 
higher  concentration  (i.e.  a  local  maximum  concentration),  but  to  send  the  sensors  to  areas  of  higher  state 
estimation  error.  An  overview  of  the  proposed  scheme  is  presented  in  Figure  1. 

2  Physical  model 

The  physical  model  employed  for  the  estimation  theory  developed  under  FA95  50-09- 1-0469  is  briefly 
described.  The  gaseous  release  under  investigation  is  assumed  to  occur  on  a  2D  plane  given  by  D  = 
[0 ,Lx\  x  [0,  Ly]  C  M2  which  is  assumed  to  be  oriented  parallel  to  the  surface  of  the  earth  [1,  2,  3].  A 
moving  point  source  S(t,  X,  Y)  =  b(X.  Y)u(t )  is  considered,  having  a  release  rate  u(t)  and  a  spatial  dis¬ 
tribution  b(X,  Y)  with  the  spatial  variables  (X,  Y)  G  Q  and  represented  as 

b(X,Y)  =  S(X-Xc(t))6(Y-Yc(t)),  (X,Y),(Xc(t),Yc(t))€n,  t  6  M+  (1) 

where  the  time  varying  centroid  of  the  source  is  denoted  by  9c{t)  =  (Xc(t).  Yr(t))  G  Q  [4,  5].  It  is  as¬ 
sumed  that  the  dispersion  of  the  gas  is  in  the  turbulent  diffusion  regime,  that  there  are  no  chemical  reactions, 
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the  ambient  speed  (wind)  is  constant,  and  that  there  is  no  impact  on  the  incompressible  background.  Un¬ 
der  these  conditions  the  dispersion  can  be  described  by  the  2D  advection-diffusion  equation  for  the  mean 
concentration  c(t,  X ,  Y)  given  by  [1,  2,  3] 


ft  +  Wxm  +  Wyw  =  w  {Kxxw) +  w  + 

c  (0,  X,  Y )  =  c0  (X,  Y) ,  {X,  Y)£to,  t  E  M+ 

where  Wx,  Wy  are  the  mean  (wind)  speed  of  the  background,  and  Kx ,  Ky  are  the  eddy  diffusivities  in 
the  X  and  Y  directions  respectively.  The  resulting  equation  with  the  Dirichlet  boundary  conditions  and 
unknown  initial  conditions  is  expressed  as 


dc  TTr  dc  Tjr  dc  d2c  d2c  .  . 

m  +  Wxox  +  Wrm  =  KxxdTi  +  Krrm+s(t'Sc)' 

c  (0,  X,  Y)  =  c0  (X,  Y) ,  c  ( t ,  0,  0)  =  c  (f,  Lx,  0)  =  c  (t,  0,  Ly)  =  c  (t,  Lx,Ly )  =  0. 


(2) 


3  Plant  data  generation 


In  the  absence  of  measurements  taken  with  the  SAV,  the  sensor  data  (i.e.  the  plant)  must  be  generated  using 
(2).  The  numerical  solution  of  (2)  is  obtained  with  a  finite- volume  method  (FVM)  [6].  The  basic  steps  of  this 
method  are  summarized  following  the  implementation  outlined  in  Gatsonis  et  al.  [7] .  The  computational 
domain  Q  is  discretized  with  a  set  of  (Nx  + 1)  x  (Ny  +  1)  grid  points  defining  the  centers  of  X  =  Nx  x  Ny 
rectangular  finite  volumes.  The  PDE  (2)  is  written  in  conservative  form  as 


dc  d(cWx )  d(cWy)  d  f  dc  \  d  f  dc  \ 
dt+  dX  +  OY  dX  \  xdX )  ~  dy  \  ydYj 


S(t,9c). 


(3) 


This  flux- form  is  integrated  over  a  fixed  finite  volume  ilc  with  a  surface  area  A  =  An  as 


d_ 

dt 


F-n  dS 


SdV, 


where  the  total  flux  term  is  expressed  as  the  sum  of  the  convective  and  diffusive  terms 


(4) 


F  =  (/£  +  /£)x  +  (/£  +  /£)Y, 
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Figure  2:  Graphical  representation  of  the  state  error  for  guidance  scheme. 


with 

(Q/i  r)n 

fx=cWX,  fy  =  cWY  f°  =  -Kxx  —  ,  fV  =  -KyY  —  . 

The  integration  (4)  over  a  cell  volume  Qc,ij  with  vertex  (ij)  shown  in  Figure  2  results  in  the  semi-discrete 
equation 

dcij/dt  =  FT),,,,  £  (Fg.A g  +  +  Fg.Ag  +  Fgfc.Agfc)  +  Sid.  (5) 

sides 

The  flux  evaluation  requires  in  general  the  local  normal  to  a  surface.  For  example,  on  the  east  surface  of  the 
volume  ( ij )  the  normal  is 

nfj  =  nEx,i3±  +  nYijY, 


and  the  flux  is 


FgAg  =  (f§E  +  >yAg  +  (gEE  +  g^:)n^A 


ij- 


For  the  structured  grid  with  rectangular  finite  volumes  considered  here,  the  local  normal  is  positive  away 
from  each  finite  volume  surface 


n  g  =  X,  n#  =  -X,  ng  =  Y,  ng  =  -Y, 

^Fjk-Afjk  =  uFjk  +  fijk)Afjk- 
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The  convective  flux  f-jE  at  the  east  boundary  is  calculated  as  the  average  of  cell  centered  fluxes  for  adjacent 
cells  as 


fCE  _  1  /  rCE  ,  fCE\ 

Jij  2  '  J  ij  )  ' 

The  diffusive  flux  for  the  east  boundary  is  calculated  with  a  central  difference  approach  at  the  cell  interface 
as 


j?DE 

Jij 


Kxx 


dc 

dX 


E 


Kxx  | 


Q+1J  Qj 

'iXi+u-Xij 


The  flux  calculation  for  the  remaining  three  sides  follows  similarly.  The  cell  Peclet  number  is  a  measure 
of  the  relative  strengths  of  the  convective  and  diffusive  parts  of  the  flow  and  is  defined  as  Pe  =  W/K.  At 
larger  Peclet  numbers  (\Pe\  >  2),  the  convective  fluxes  are  evaluated  using  upwinding.  For  Wx  >  0, 


fCE  _  fCE 
J ij  J  ij 


cwx Uj,  .ii;w 


fCW 
J i—l,j 


cWx  | 


j-ij' 


For  Wx  <  0, 


fCE  _  fCE  _  „u/  |  fCW  _  fCW  _  | 

Jij  ~  Ji+lJ  ~  °W X  I  j+1  j  i  Jij  ~  Jij  ~  CWx | 


The  semi-discrete  equation  (5)  can  be  written  as 


dcij/dt  —  [cj+ijdjj  T  Cj— \ ,j(Jjj  Cij+i^ij  T  P.j— l  Ay  T  T  Sj,^ 


(6) 


The  coefficients  are  given  by 


w 

a-  =  max 


F\v,  (  Dw  +  )  ,0 


aE:  =  max 


-FeADe-^  )  ,0 


and  D\y,  De,  F\y  and  FE  are  defined  as 

Dw  =  j—^—Aw,  De  =  VE  A E,  Fw  =  (cWx)wAw,  Fe  =  (cEx)eAe. 

0A.wp  0A.pp 

where  5X\yp  is  the  distance  between  the  volume  (ij)  and  the  west  volume.  Application  of  (6)  to  all  cell 
centers,  leads  to  a  system  of  ( Nx  x  Ny )  =  N  semi-discrete  ODEs. 

Introducing  x  =  c  =  {ci, . . . ,  cE }  as  the  finite  dimensional  representation  of  the  state  vector,  the 
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indexing  for  mapping  from  (ij )  to  n  is  expressed  as 


Cn  =  Cij,  n  =  U  -  l)Nx  +  i,  j  =  1,  Ny,  i  =  1,  Nx- 


The  system  of  semidiscrete  ODEs  can  then  be  written  in  the  state-space  form  as 


x  =  Tlx  +  Bu. 


(V) 


The  source  is  represented  by  the  release  rate  u  (t)  and  the  spatial  distribution  associated  with  the  source 
location  (Xc,  Yr),  (1).  For  a  point  source  release  located  within  volume  (ij),  the  vector  B  is  all  zeros  with  a  1 
located  at  B  ( n ).  The  (N  x  N)  state  matrix  A  is  associated  with  the  finite  dimensional  approximation  of  the 
advection  diffusion  equation.  The  operator  A  ( t )  is  constructed  with  the  a  coefficients.  The  ap  coefficient 
located  at  A  ( n ,  n)  is  the  entry  for  the  nth  finite  volume  and  thus  the  entry  for  the  corresponding  state 
x„  .  The  remaining  coefficient  are  located  at  A  (n,  m),  where  m  is  the  index  location  of  the  corresponding 
neighboring  node. 

The  system  of  semi-discrete  ODE's  as  in  equation  (5)  is  integrated  using  the  four-step  Runge-Kutta 
scheme  expressed  as 

CM  =  c(n)  _  amA tR(AV)c(m-i)^  m  =  1,2,  3, 4,  c(n+1)  =  c(4) 

where  is  the  concentration  at  time  level  n  and  a\  =  0.25,  «2  =  0.333,  a 3  =  0.5,  0:4  =  1.0. 

3.1  Sensor  measurements  from  numerically  generated  process  model 


The  sensor  is  assumed  to  have  measurement  data  on  the  concentration  and  gradient  at  its  current  spatial 
location.  Measurement  data  is  taken  from  the  simulation  of  the  plant.  The  state  measurement  at  the  spatial 
location,  c(X,  Y),  is  calculated  as  the  weighted  average  based  on  the  inverse  distance  to  the  4  nearest  volume 
centers 


c(X,Y)  = 


CNE  *  dNlE  +  CSE  *  dslE  +  csw  *  dgly  +  CNW  *  dNlw 

dxE  +  dSE  +  dsw  +  dNW 

dc  dc 


The  structured  discretization  of  the  grid  allows  the  local  gradient  (-§x,  to  be  calculated  with  the  second 
order  central  differencing  scheme  for  nonuniform  spacing  [6,  8].  This  approach  calculates  the  gradient  based 
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on  two  computational  points  in  each  direction  and  the  distance  to  the  points. 


4  Model-based  State  Estimation  and  Sensor  Guidance 

4.1  Output  measurement  model 

Numerous  sensors  are  available  for  the  measurement  of  trace  species  in  the  atmosphere.  They  are  available 
in  several  types  with  many  different  operating  characteristics  including  size,  sensitivity,  dynamic  range, 
reliability,  and  response  time  [9,  10,  11,  12].  Here,  a  specific  sensor  and  contaminant  pair  have  not  been 
chosen,  so  a  generalized  sensor  model  is  constructed  based  on  realistic  sensor  parameters.  It  is  assumed  that 
the  sensor  is  selective  and  able  to  distinguish  the  contaminant  from  the  other  constituents  in  the  atmosphere. 

A  mobile  agent  equipped  with  a  sensor  is  able  to  provide  concentration  information  c(f ,  X,  Y)  at  the 
spatial  point  6S  =  (Xs,  Ys)  G  Q.  Similar  to  the  source  model  in  (1),  the  spatial  distribution  of  the  sensor 
is  also  modeled  by  a  spatial  Dirac  delta  function.  It  is  assumed  that  there  is  no  noise  and  the  measurement 
device  provides  exact  readings  of  the  local  concentration 

y(t;9s)  =  c(t,Xs,Ys)=  6(X  -  XS)5(Y  -  Ya)c(t,  X,Y)  dX  dY 

Jo  Jo 

With  the  sensor  mounted  on  a  mobile  agent  (S  AV),  it  is  able  to  provide  measurement  data  at  various  locations 
within  the  spatial  domain  9  as  a  function  of  time.  The  time  dependent  sensor  measurement  taken  at  the 
centroid  0s(t)  of  the  sensor  is  then 

y(t-es(t))=  [Lx  [LY  8{X  -  Xs(t))6(Y  -Ys(t))c(t,X,Y)dX  dY.  (8) 

Jo  Jo 

which  is  essentially  the  concentration  c(t,  X,  Y)  at  the  current  position  0s{t)  of  the  SAV. 

4.2  Abstract  formulation  of  advection  diffusion  PDE 

The  advection  diffusion  PDE  (2)  can  be  viewed  as  an  evolution  equation  in  an  appropriate  Hilbert  space. 
Such  an  abstract  formulation  is  conducive  to  both  the  ensuing  stability  analysis  for  the  resulting  state  esti¬ 
mator  and  its  finite  dimensional  approximation. 

The  state  space  is  taken  to  be  A  =  L2(Q)  equipped  with  inner  product  denoted  by  (•,  ■)  and  norm  by 


|  •  |.  Additionally,  we  consider  V  =  Hq  (Q)  be  the  Sobolev  space  with  V  dense  in  X.  Such  a  consideration  is 
necessitated  as  the  output  and  input  operators  are  defined  in  V  and  its  dual  V*  [13].  The  concentration  state 
is  an  element  of  the  Hilbert  space  x(t)  =  c(t,  ■,  ■)  in  X  over  [0,  T],  The  PDE  (2)  can  then  be  written  as  an 
evolution  equation  [14,  15] 


x(t)  =  Ax(t )  +  B(9c(t))u(t),  x(0)  =  xq  €  X 

y(t;0s(t ))  =  C(9s(t))x(t) 


(9) 


where  A  €  C(V.  V*)  is  the  elliptic  operator  associated  with  the  advection  diffusion  PDE  (2),  B{9c(t))  is 
the  location-parameterized  input  operator  associated  with  the  source  spatial  distribution  b(X,  Y)  in  (1),  and 
CUhit)  )  is  the  output  operator  associated  with  the  time  dependent  sensor  measurement  (8).  For  the  problem 
to  be  well  posed,  one  needs  to  impose  B(9c(-))f(-)  €  L2( 0,  T;  V*),  see  [16,  17,  15].  Equation  (7),  as  used 
for  plant  generation,  constitutes  a  finite  dimensional  representation  of  the  evolution  equation  (9). 


4.3  Model-based  State  estimation 

The  estimator  developed  in  [18]  is  modified  to  account  for  vehicle  dynamics  and  is  applied  in  this  work 
to  estimate  the  concentration  state  over  the  entire  spatial  domain.  In  summary,  a  Luenberger  observer  is 
implemented  with  the  filter  gain  taken  to  be  a  constant  multiple  of  the  dual  of  the  output  operator  (collocated) 
and  given  by 

x{t)  =  ^A-jC*(9s(t))C(9s(t))'jx(t)  +jC*(9s(t))y(t;9s(t)),  x(0)  /  x(0),  7  >  0,  (10) 

Equation  (10)  constitutes  half  of  the  integrated  estimation  and  guidance  equations.  It  is  accompanied  by 
another  equation  that  provides  the  time-variation  of  the  sensor  centroid  9s(t),  which  would  dictate  the  spatial 
repositioning  of  the  SAV  within  the  spatial  domain  1>.  This  spatial  repositioning  is  given  implicitly  in  terms 
of  the  control  inputs  (torques)  to  the  aerial  vehicles  that  carries  the  sensor. 

Sensor  guidance  is  dependent  on  the  state  estimation  error,  that  is  e(t)  =  x(t)  —  x(t).  The  spatially 
distributed  state  error  e(t,  X,  Y)  =  x(t,  X,  Y)  —  x(t,  X,  Y)  may  be  denoted  simply  as  e(t)  when  considered 
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as  an  element  of  the  Hilbert  space  X.  Using  (9)  and  (10),  the  evolution  of  the  state  error  is 


e(i)  =  A(9s(t))e(t)  +  B(9c(t))u(t), 
e(0)  =  e0  €  X. 


wher  eAd(0s(t))  =  (A-'y  C*(9s(t))C(9s(t))j. 


(11) 


4.4  Sensing  device  model 

The  physical  construction  of  sensors  introduces  inherent  limitations  in  the  sensor  operation.  Ram  et  al  [9] 
provide  a  good  discussion  on  sensor  construction  and  design.  Many  trace  contaminant  sensors  do  not  take 
continuous  measurements,  but  instead  have  a  recovery  time  or  transient  response  time  between  sensor  read¬ 
ings  [10,  19,  20].  Transient  effects  limit  the  time  between  accurate  sensor  readings.  A  sensor  that  provides 
accurate  measurements  at  a  very  fast  sample  rate  will  clearly  provide  more  information  than  one  that  sam¬ 
ples  slower.  This  research  effort  is  interested  in  real  time  state  estimation,  so  faster  measurement  frequencies 
are  desirable.  Optical  based  sensors  such  as  chemiluminescence  detectors  and  absorption  spectroscopy  are 
able  to  sample  every  1-5  seconds  [10].  For  the  numerical  experiment  considered  in  this  research  effort,  time 
between  samples  is  taken  to  be  2  seconds. 

Sensor  dead-zone  and  threshold  saturation  are  also  implemented  in  the  model.  All  trace  contaminant 
sensors  have  maximum  and  minimum  sensing  thresholds  that  determine  the  limits  of  the  sensing  ability. 
When  the  state  is  below  a  minimum  threshold  cmin,  the  sensing  device  does  not  detect  an  elevated  measure¬ 
ment.  This  is  referred  to  as  a  dead-zone.  When  the  state  is  above  the  maximum  cmax,  the  output  from  the 
sensor  is  the  saturated  value.  These  limits  are  expressed  below 

0  <C  cmin 

c(Xs,Ys ) 

Cmin  <  c(Xs,Ys)  <  Cmax  (12) 

Cmax  c(Xs,Ys)  > 

Cmax 

When  taking  sensor  measurements  on  the  environment,  the  time  rate  of  change  of  the  sensor  depends 
on  the  process  under  investigation  (the  atmospheric  advection  diffusion)  as  well  as  the  motion  of  the  sensor 


cS(Xs,Y, )  = 
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and  the  spatial  gradient  of  the  concentration  [21]  as  shown  below 


dc  ( dc\  dX  (  dc  \  dY  (  dc  \ 

dt  \^J  x,y  +  dt  \dxJY,t+  dt  \dY )  x,t 

where  ^  is  the  velocity  of  the  sensor  in  the  X  direction  and  is  the  spatial  gradient  of  the  contaminant 
in  the  X  direction.  This  effort  simplified  sensor  measurements  and  assumed  sensor  readings  at  the  current 
location  of  the  sensing  agent  were  directly  available. 

4.5  SAV  dynamics 

The  movement  of  the  sensor  throughout  the  domain  is  accomplished  with  an  SAV.  For  simplicity,  the  sensor 
location  is  taken  to  be  at  the  barycenter  of  the  SAV.  The  equations  of  motion  describing  the  mobile  agent 
position  will  then  be  valid  for  the  position  of  the  sensor  centroid. 

Depending  on  the  focus  of  the  work,  mobile  agents  may  be  modeled  with  great  detail  [22,  23]  or  very 
simply  [24,  25,  26,  27].  In  this  effort,  the  mobile  agent  is  modeled  as  a  fixed  wing  unmanned  aerial  vehicle 
with  basic  autopilot  controls.  With  knowledge  of  its  own  state,  the  lower  level  controllers  account  for 
deviations  in  the  trajectory  of  the  aircraft  caused  by  wind,  coupled  control  surfaces,  and  other  disturbances. 
The  SAV  is  assumed  to  be  rigid  and  symmetric  with  a  collocated  center  of  mass  and  center  of  gravity,  which 
allows  the  equation  of  motion  to  be  presented  in  a  more  compact  form.  The  mass  and  moment  of  inertia  are 
assumed  constant  throughout  the  simulation.  For  simplicity,  many  works  focus  on  the  kinematic  motion  and 
neglect  the  dynamics.  However,,  in  order  to  provide  a  more  accurate  representation  of  the  SAV  motion,  this 
research  effort  considered  the  dynamic  equations  of  motion. 

The  SAV  motion  is  constrained  due  to  the  physical  limitations  of  the  aircraft.  The  agent’s  speed,  v(t), 
and  turning  rate,  are  constrained  within  maximum  and  minimum  values  [27,  28,  29]. 

0  <  Vmin  <V<  Vmax,  ~U^max  —  max ' 

The  thrust  and  turning  force  are  also  constrained  as 


~Tlr 


<  n  <  Tir 


~T, 


T  T 

a  max  —  ‘a  •  a  max  • 


For  an  aircraft  equipped  with  a  lower  level  autopilot  [29,  23],  the  pose  of  the  sensor  may  be  described 
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V 


Figure  3:  Schematic  of  the  SAV  coordinate  system  in  2D. 


by  the  kinematic  equations  [27]  given  by 


X(t)  =  v(t)cos(i/;(t)),  Y(t)  =  v(t)  sin('i/’(i)),  ip(t)  =  (13) 


where  the  pose  is  represented  by  the  cartesian  position  (X,  Y )  and  the  SAV’s  heading  angle  A.  The  motion 
is  determined  by  the  SAV’s  speed  and  angular  turning  rate  (v(t),cu^(t))  as  shown  in  Figure  3.  Equation  (13) 
can  be  expressed  in  matrix  form 

q(f)  =  S  (t)v(f)  (14) 


where  q(i) 


X(t),Y(t),ip(t)  ,  v(f) 


T 


and  the  transformation  matrix  S (t)  is 


m 


cos  (ip)  0 
sin(^)  0 

0  1 


(15) 


These  equations  are  the  same  for  those  of  a  differential  drive  robot  [29],  with  the  additional  constraint  that 
the  fixed  wing  aircraft  must  hold  a  minimum  forward  velocity  to  maintain  lift. 

The  2D  dynamic  equations  of  motion  for  steady  state  flight  given  in  inertial  coordinates  [30,  31,  24,  25, 


12 


23]  are  given  by 


X 

COS  (tjj) 

0 

—M'fvsin'f 

n 

Y 

sin  (ip) 

0 

Ta 

+ 

Mfvcos'f 

$ 

0 

1 

0 

where  77  is  the  thrust  and  ra  is  the  angular  torque  applied  by  the  control  surfaces  to  the  center  of  mass. 
Since  the  equations  are  for  steady  flight,  there  is  no  explicit  drag  component.  With  77  and  ra  set  to  zero, 
the  aircraft  maintains  constant  speed  and  angular  turning  rate.  The  mass  matrix  A4  is  expressed  as  Ai  = 
diag{M,  M,  /},  where  M  is  the  mass  and  /  the  moment  of  inertia  of  the  SAV.  Equation  (16)  can  be  rewritten 
in  matrix  form  as 

Afq  =  B\t  —  ArX.  (17) 

Substituting  equation  (14)  into  equation  (17),  along  with  q  =  Sv  +  Sv,  result  in 

M( Sv)  =  Bxr  -  AtX  =>  M Sv  +  A4Sv  =  Bit  -  AT\. 

Premultiplying  by  ST  yields 


StMSv  +  StMSv  =  StB1t-StAt\.  (18) 

The  equation  is  simplified  with  the  introduction  of  the  transformed  mass  matrix  =  S7  A4S  =  diag{M,  I } 

The  term  SrA4 1 S  is  simplified  to  SrA4iS  =  02xi-  The  input  control  matrix  B>  is  also  simplified  to 
B'2  =  STB\  =  12x2-  As  a  consequence  of  the  premultiplication  by  ST,  the  constraint  matrix  is  eliminated 
from  (17)  since  S7  AA  =  02xi-  For  the  inputs  r  =  [77.  r„],  the  dynamic  equation  of  motion  from  (18)  can 
then  be  written  as 

A4iv  =  B2T.  (19) 

Remark  1  When  the  SAV  is  assumed  massless  and  inertialess,  the  kinematic  equations  (14)  can  be  used  to 
describe  the  SAV  motion  in  2D.  If  it  is  further  assumed  that  its  motion  is  holonomic  and  not  constrained  by 
(14),  then  it  is  able  to  implement  motion  with  any  desired  direction  and  speed.  In  fact,  this  was  considered 
in  [18]. 
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4.6  SAY  guidance 


The  Lyapunov  functional  used  in  [18,  32,  33]  is  modified  to  include  the  kinetic  energy  of  the  SAV 

V  =  -(e(t),Ad(Oa(t))e(t))+KE. 

The  kinetic  energy  of  SAV  is  related  to  its  mass  and  velocity  as  KE  =  jqT.Adq.  Transforming  from 
cartesian  velocities  to  body  velocities  with  equation  (14)  yields  KE  =  ^vTSi  A4Sv  =  M.\v.  For  an 

aircraft  with  its  motion  constrained  in  a  2D  plane  parallel  to  the  ground,  the  gravitational  potential  energy 
term  is  equal  to  zero.  The  Lyapunov  functional  is  re-written 


V  —  —  (e,^4  ci(Os)e)  + 


1 

2 


vtA4iv 


(20) 


Taking  the  time  derivative  of  the  Lyapunov  functional  along  the  trajectories  of  the  state  error  (11)  and  the 
complete  dynamics  of  the  craft  (14),  (19)  yields 


V 


—  |.Acz(0s)e|  +  EEx A"  +  EEyY  I  +  V  A4iV 


—  \Ad{0s)e\  + 


EEx  0 


q  +  v  A4iv. 


Substituting  equations  (14)  and  (19)  simplifies  to 

V  =  -\Ad{es)e\2  + 
=  ~  |“4ci(@s)e|2  + 


eex  0 

eex  ££y  0 


Sv  +  vrM  i  v 
Sv  +  vT(32r; 


(21) 


where  e  denotes  the  state  estimation  eiTor  at  the  current  sensor  location  and  Ex ,  £y  denote  the  spatial  gra¬ 
dients  of  the  estimation  error  at  the  sensor  location.  The  control  law  can  now  be  developed  to  ensure  that 
the  Lyapunov  derivative  is  negative  semi-definite  and  the  guidance  scheme  drives  the  sensor  to  areas  of 
higher  state  error.  Focusing  on  the  part  of  equation  (21)  to  be  made  negative  definite  and  recognizing  that 
vT  Bot  =  (vt£>2t)T  gives 

eex  ££y  0 


Sv  +  {B2t)t  v  <  0 
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Solving  for  the  input  torque  yields  the  following  control  guidance  law 


lc  T 


T  =  -£>2  S 


EEx  ££y  0 


To  better  understand  the  control  law  (22),  the  matrices  are  expanded  as 


(22) 


- 

Tl 

1 

II 

1 

1 _ 

1 

COS  f 

0 


sin  -0 

0 


0 

1 


EEy 

0 


— EEx  COS  f  —  EEy  sin  if) 

0 


(23) 


The  above  relates  the  SAV  motion  control,  via  the  torque  r,  to  the  performance  of  the  estimator,  in  other 
words  the  motion  of  SAV  is  dictated  by  the  performance  of  the  estimator  and  is  given  explicitly  in  terms  of 
the  output  estimator  error  eU)  and  the  spatial  gradients  Exit),  Ey(t  )  of  the  estimation  error  at  the  current 
sensor  location. 

The  control  torque  can  further  be  modified  to  make  V  =  —  \A,:i(0fe\2  —  vv  /v  v  “more”  negative  by 
including  velocity  terms 


ti  =  —ci  (ee x  cos  f  +  EEy  sin  f)  —  k\v  —  ki2ip,  ci  >0,  k\  >  0,  k\2  >  0 

(24) 

ra  =  —k\2V  -  k2lp,  k2  >  0,  ki2  >  0. 

The  constants  c\,k\,  k\2  are  user-defined  guidance  gains,  chosen  to  achieve  desired  sensor  performance.  In 
this  research  effort,  ci  was  chosen  so  that  the  input  torque  is  of  the  same  order  of  magnitude  as  the  agent 
physical  limitations  (i.e.  the  input  is  not  always  saturated).  The  —  k\v  —  k\2’f  term  was  not  required  for  this 
instance. 

The  Lyapunov-based  guidance  scheme  (22)  only  provides  motion  control  inputs  for  the  thrust  77.  While 
the  extension  in  (24)  could  possibly  include  a  velocity  term  in  the  rotational  dynamics  via  the  additional 
entries  of  the  positive  semi-definite  matrix  K 


If  =  Ta  =  —  k\2V  —  k2f,  k2  >  0,  k\2  >  0 

it  is  nonetheless  unable  to  address  the  rotational  motion  in  an  effective  manner.  A  supplementary  controller 
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was  used  whereby  the  angular  torque  input  was  given  by 


Ta(t)  =  I^d{t)  -  kd  -  kP  (VM  -  V,d(i))  (25) 

The  current  heading  angle  was  taken  to  be  ipa(t).  The  desired  heading  angle  ipd(t)  was  defined  as  the  angle 
created  by  the  two  cartesian  torques 

i/jd  =  arctan  (  —  ) 

VxJ 

The  above  choice  stems  from  the  work  in  [18]  which  considered  a  massless  and  inertialess  sensing  agent. 
From  equation  (25),  ipd,  the  time  derivative  of  the  desired  angle  was  calculated  numerically. 

4.7  SAV  deployment  model 

Prior  to  initiating  the  tracking  scheme,  the  SAV  assumes  no  source  is  present  in  the  domain.  With  knowledge 
of  the  current  wind  profile,  the  sensor  is  strategically  placed  downwind  in  the  domain.  As  the  wind  blows, 
it  will  advect  contaminant  towards  the  sensor.  Although  patrolling  downwind  of  the  source  will  increase 
the  probability  of  detecting  a  source  in  the  domain,  the  search  strategy  is  not  assumed  optimal.  The  sensor 
travels  in  a  circular  path  downwind  continuously  taking  readings  until  it  detects  an  elevated  concentration, 
(12).  The  sensing  agent  then  stops  patrolling  and  begins  searching  for  the  source  via  (10),  (22),  (25).  At  this 
point,  the  detection  scheme  begins  estimating  the  process  state  and  providing  the  mobile  sensing  agent  with 
the  appropriate  command  signals  from  the  guidance  scheme. 

5  Finite  dimensional  approximation  and  sensor-based  grid  adaptation 

The  estimator  is  approximated  as  a  finite  dimensional  system  similar  to  the  approach  used  with  the  plant 
in  (7).  A  conservative  form  of  the  advection  diffusion  equation  is  used  and  integrated  with  the  four  stage 
Runge-Kutta  scheme.  Since  it  is  desired  to  calculate  the  estimator  in  real  time,  the  dimension  of  its  finite 
dimensional  approximation  is  reduced  compared  to  that  of  the  plant.  To  maintain  a  desired  level  of  accuracy 
in  regions  of  interest,  while  keeping  the  computational  requirements,  grid  adaptation  is  implemented.  The 
implementation  of  a  reduced  dimensional  estimator  also  avoids  the  inverse  crime  problem  that  is  associated 
with  a  numerical  simulation  and  inversion  that  are  of  the  same  discretization.  An  equivalent  discretization 
in  the  forward  simulation  and  estimated  state  would  provide  results  that  are  overly  optimistic  [34]. 
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Figure  4:  Graphical  representation  of  the  state  error  for  guidance  scheme. 

The  limited  computation  power  and  desire  for  a  real  time  solution  suggested  the  use  of  grid  adaptation 
techniques  [35,  36].  For  this  research  effort,  a  combination  of  a  stretched  grid  and  switched  system  were 
implemented.  The  estimation  scheme  discretized  the  domain  with  9  stretched  grids.  Each  of  the  stretched 
grids  consisted  of  a  predefined  number  of  volumes  in  each  direction.  This  created  a  system  that  was  low 
enough  in  dimension  to  be  solved  in  real  time.  Each  stretched  grid  had  an  area  of  relatively  high  resolution 
that  was  meant  to  focus  on  the  area  of  interest.  The  rest  of  the  domain  consisted  of  a  coarser  discretization 
to  keep  computational  requirements  low.  Autonomous  state  dependent  switching  [37]  control  was  applied. 
When  the  area  of  interest  changed,  the  grid  was  switched  to  ensure  the  location  of  interest  always  had  a 
higher  degree  of  discretization  than  the  rest  of  the  domain. 

This  hybrid  dynamical  system  couples  the  estimation  scheme  with  the  computational  scheme,  using 
one  to  enhance  the  other.  Numerically,  this  switching  changes  several  of  the  matrices  used  in  the  state 
space  formulation.  At  each  switching  instance,  the  state  matrix  A  had  to  be  recalculated.  The  data  from 
the  old  stretched  grid  was  moved  to  the  new  one  by  means  of  a  prolongation  and  restriction  operation. 
A  prolongation  operator  transferred  information  from  the  old  grid  to  a  new,  higher  dimensional  one  with 
higher  resolution.  A  restriction  operator  then  transferred  the  data  to  the  new  grid.  Nine  restriction  and 
prolongation  operators  were  created  a  priori  in  order  to  handle  all  switching  cases  [35].  This  helped  to 
minimize  computations  during  the  state  estimation.  The  output  C  matrix  changed  just  slightly  due  to  the 
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Ai  A5  A8 

Figure  5:  Example  of  grid  adaptation  with  high  resolution  area  highlighted. 


spatial  location  of  the  sensor,  but  since  it  was  created  at  each  time  step,  information  did  not  have  to  be 
transferred  from  one  to  another  when  the  switch  occurred.  Below,  Algorithm  1  outlines  the  grid  adaptation 
logic.  Figure  5  demonstrates  the  grid  adapting  from  Grid  1,  to  Grid  5,  then  to  Grid  8  as  the  sensor  position 
moves  between  the  locations  noted.  The  set  of  available  switched  grids  resulted  in  a  family  of  matrices 


Algorithm  1  State  dependent  grid  adaptation 

1:  estimate  (X,  Y,  ij))c 

2:  calculate  nearest  switched  grid 

3:  if  nearest  grid  /  current  grid  then 

4:  set  new  grid  to  nearest  grid 

5:  switch  state  matrix,  A 

6:  prolongate  concentration  data  to  general  grid 

7:  restrict  state  information  to  new  grid 

8:  end  if 


{Ai,  i  €  1}  based  on  the  index  set  X.  In  this  research  effort,  the  switching  signal  (which  was  based  on  the 
sensor  location)  was  a  piecewise  constant  function  in  time  represented  by  a  :  [0,  oo)  — »•  X.  Consider  the 
family  (Sp)p&z  of  linear  systems  that  are  continuous  in  time.  For  each  p  €  X,  the  estimator  is  given  by 

xn{t)  =  (Ap  -  7 Cp  (es(t))Cp{9s(t)))  xn{t)  +  7 Cp  (0s(t))y(t;  9s{t ))  (26) 

The  current  sensor  position  9s{t)  dictates  the  switching  signal  and  subsequently  the  choice  of  the  matrices 
Ap,  Cp.  At  the  same  time,  the  state  estimator  along  with  the  guidance  scheme  (10),  (22)  provide  the  spatial 
repositioning  of  the  SAY. 
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Figure  6:  Simulation  Flow  Chart. 


6  Implementation  pseudo  code  and  flowchart 


Due  to  the  lack  of  experimental  data,  simulations  for  this  work  were  carried  out  in  two  parts.  The  first  part 
consisted  of  generating  acceptable  sensor  data  that  the  estimator  could  access  in  place  of  real  measurements 
as  described  in  Section  3.  A  source  in  a  2D  domain  was  simulated  on  a  high  dimensional  system  with  the 
concentration  profile  saved  to  file  at  every  time  step.  The  high  dimensional  approximation  was  computation¬ 
ally  expensive  and  was  therefore  done  a  priori  and  stored  to  file.  In  practice,  this  step  would  not  be  required 
and  measurement  data  would  be  taken  directly.  The  pseudocode  for  the  generation  of  experimental  data  can 
be  seen  in  Algorithm  2.  The  actual  estimation  was  done  in  the  second  part  of  the  simulation.  As  the  sensing 


Algorithm  2  Generating  of  Sensor  Data 

1:  read  simulation  parameters 
2:  discretize  high  dimensional  uniform  grid 
3:  for  t  =  dt  to  t final  do 
4:  calculate  source  location  ( X ,  Y)c 

5:  RK4  integration  of  forward  state  x 

6:  apply  boundary  conditions 

7:  end  for 

8:  output  state  at  each  time  step  to  file 


aerial  vehicle  moved  around  the  spatial  domain,  it  was  continuously  taking  measurements,  which  consisted 
of  reading  small  portions  of  the  stored  data  from  the  files  created  in  the  generation  of  source  data.  Algorithm 
3  details  the  estimation  process  and  the  entire  simulation  is  outlined  in  Figure  6. 
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Algorithm  3  State  Estimation  Scheme 
Require:  Output  data  files  from  forward  problem. 
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end  if 

RK4  integration  of  x 
apply  BCs 
else 

continue  patroling 

end  if 

calculate  new  ( X ,  Y)s 
switch  grid 

end  for 
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7  Summary  of  Simulation  Results 

Several  simulations  have  been  performed  with  the  developed  guidance  scheme.  For  consistency,  the  domain 
size  in  each  case  was  taken  to  be  4  km  x  4  km.  The  eddy  diffusivity  was  assumed  constant  in  space  and  time 
as  20m2/s.  The  wind  was  also  assumed  constant  in  space  and  time,  blowing  at  5m/s  to  the  East  and  5m/s 
to  the  North.  In  all  cases,  the  sensor  started  in  the  downwind  area  of  the  domain  in  a  patrolling  behavior. 
Simulations  are  conducted  on  a  5  node  Linux  cluster  running  Red  Hat  3.4.6.  The  serial  code  was  implanted 
on  one  of  the  nodes  with  a  Quad  Core  Intel  Xeon  processor  running  at  2.33GHz  and  16GB  of  RAM.  The 
code  was  compiled  with  Intel’s  Fortran  compiler  IFORT  version  9.0.  Data  visualization  was  performed  with 
Tecplot  360  software  version  12  running  on  a  Windows  7  workstation. 

The  numerical  discretization  of  the  forward  problem  (plant)  in  each  case  was  chosen  to  be  a  90  x  90 
structured,  uniform  grids.  The  estimator  discretization  was  significantly  coarser  at  just  30  x  30  computational 
points.  These  discretizations  ensured  high  fidelity  data  was  available  from  the  forward  problem,  while  still 
allowing  the  estimator  to  be  calculated  in  real  time.  A  computational  time  step  of  0.1s  was  chosen  for  all 
simulations,  which  was  significantly  below  the  value  dictated  by  stability  requirements. 
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Figure  7:  Trajectory  of  the  sensor  (green)  and  with  a  stationary  source  (red). 

7.1  Stationary  Source 

A  stationary  source  was  simulated  with  a  constant  release  rate  in  the  center  of  the  domain.  Such  a  stationary 
source  would  be  expected  for  a  crash  or  leak  situation.  The  entire  500s  simulation  took  187s  to  calculate. 
Figure  7  shows  the  trajectory  of  the  sensor  as  it  travels  towards  the  source.  At  approximately  110s,  the 
sensor  detected  a  nonzero  concentration  and  begun  the  estimation  process.  As  shown  in  Figure  8,  the  sensor 
got  very  close  to  the  source  in  approximately  180s.  From  the  plot  of  the  norm  of  the  concentration  error 
given  in  Figure  9,  it  can  be  seen  that  the  norm  is  increasing  until  the  sensor  starts  heading  towards  the  source. 
The  error  norm  then  quickly  falls.  The  error  then  fluctuates  around  a  non-zero  value,  as  is  expected  since 
the  source  is  stationary  and  the  sensor  can  not  stop  moving.  As  the  SAV  flies  past  the  source,  the  state  error 
at  the  location  of  the  source  increases  until  the  sensor  is  driven  back  to  the  source  location. 

7.2  Crossing  Source  Trajectory 

A  crossing  source  trajectory  was  also  considered.  The  source  traveled  directly  across  the  diagonal  of  the 
domain,  then  traveled  directly  across  the  other  diagonal.  The  source  in  this  case  had  a  maximum  velocity 
of  15m/s.  The  crossing  source  trajectory  simulated  1450s  and  took  892s  to  compute.  The  trajectory  of 
the  source  and  sensor  is  shown  in  Figure  10.  At  approximately  200s,  the  sensor  detected  the  source  and 
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Figure  8:  Distance  between  the  source  and  sensor  for  a  stationary  source. 
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Figure  9:  State  error  norm  for  a  stationary  source. 
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Figure  10:  Trajectory  of  the  sensor  (green)  and  a  crossing  source  trajectory  (red). 

begun  the  estimation  process.  Figure  1 1  shows  that  the  sensor  and  source  got  very  close  in  around  400s. 
From  the  plot  of  the  state  error  norm  in  Figure  12,  the  guidance  scheme  drives  the  state  error  towards  zero. 
Fluctuations  in  the  value  were  due  to  the  fact  that  the  sensor  was  constantly  moving. 

7.3  Circular  Source  Trajectory 

A  circular  source  trajectory  was  examined  in  the  third  simulation.  Such  a  situation  models  a  source  releasing 
material  over  a  designated  area  in  order  to  cause  damage  to  that  area.  The  circular  source  trajectory  simulated 
500s  and  took  437s  to  compute.  The  source  had  a  maximum  velocity  of  15m/s.  The  trajectory  of  the  source 
and  sensor  is  shown  in  Figure  13. 

Figure  14  shows  the  distance  between  the  source  and  sensor  as  a  function  of  time.  The  sensor  remained 
close  to  the  source  for  most  of  the  trajectory,  but  wandered  away  at  two  points.  The  first  was  around  400s. 
At  this  point,  the  sensor  went  too  far  towards  the  X-axis  and  initiated  a  spiral  outward  search  to  again  find 
the  plume.  When  it  found  the  plume,  it  traveled  downwind  where  there  was  a  high  state  estimation  error 
before  heading  towards  the  plume  upwind.  The  second  loss  occurred  around  550s.  Here,  the  sensor  traveled 
too  far  towards  the  Y  axis  and  lost  the  plume,  but  quickly  returned  toward  the  source.  The  state  error  norm 
plot  is  given  in  Figure  15.  In  both  cases  where  the  sensor  lost  the  plume,  the  norm  can  be  seen  to  increase  a 
bit.  However,  the  overall  trend  of  the  norm  is  driven  towards  zero,  as  desired. 
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Figure  1 1 :  Distance  between  the  sensor  and  a  crossing  source  trajectory. 
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Figure  12:  State  error  norm  for  a  crossing  source  trajectory. 
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Figure  13:  Trajectory  of  the  sensor  (green)  and  a  circular  source  trajectory  (red). 


Figure  14:  Distance  between  the  sensor  and  a  circular  source  trajectory. 
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Figure  15:  State  error  norm  for  a  circular  source  trajectory. 

7.4  Overlapping  Source  Trajectory 

An  overlapping  source  trajectory  was  simulated.  This  type  of  source  would  be  applicable  to  an  intruder 
doing  surveillance  over  an  area,  or  quickly  entering  an  area  to  deposit  material  and  quickly  leave.  The 
overlapping  source  trajectory  simulated  700s  and  took  226s  to  compute.  The  maximum  source  velocity  in 
this  case  was  18m/s.  The  trajectory  is  shown  in  Figure  16.  For  this  trajectory,  the  source  followed  the 
sensor  very  closely  for  the  entire  simulation.  The  sensor  detected  a  nonzero  concentration  after  120s  and 
quickly  started  following  the  source.  The  norm  of  the  state  error  for  this  trajectory  provided  some  interesting 
results.  The  norm  was  reduced  as  the  source  traveled  through  the  first  half  of  the  trajectory.  However,  as 
the  source  traveled  along  the  second  half,  it  was  traveling  downwind.  With  the  sensor  slightly  behind  the 
source  and  the  material  from  the  source  advecting  downwind,  the  state  error  norm  was  slowly  increasing, 
until  the  source  stopped  heading  directly  downwind  in  approximately  650s.  The  norm  then  again  continued 
to  decline. 


8  Conclusions 

A  model-based  estimation  scheme  of  a  mobile  gaseous  source  has  successfully  been  considered  for  a  2D 
spatial  domain.  Such  a  spatial  domain  was  assumed  to  be  parallel  to  earth’s  surface.  Using  a  mobile  sensor 
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Figure  16:  Trajectory  of  the  sensor  (green)  and  an  overlapping  source  trajectory  (red). 


Figure  17:  Distance  between  the  sensor  and  an  overlapping  source  trajectory. 
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Figure  18:  State  error  norm  for  an  overlapping  source. 

on  board  an  unmanned  aerial  vehicle  (sensing  aerial  vehicle),  the  estimation  scheme  provided  the  spatial 
repositioning  of  the  SAV  in  terms  of  control  inputs  (torques)  to  the  SAV.  Using  Lyapunov  methods,  the 
model-based  estimation  scheme  accounted  for  the  dynamic  motion  of  the  SAV  and  provided  its  guidance  in 
terms  of  the  performance  of  the  estimation  scheme;  thus  the  motion  of  the  SAV  was  solely  dictated  by  the 
performance  of  the  estimation  scheme.  Such  a  scheme  was  essentially  a  (spatial)  gradient  scheme  where  the 
SAV  was  guided  to  the  spatial  regions  of  higher  state  estimation  error  or  higher  gradients  of  the  estimation 
error. 

While  the  infinite  dimensional  model  of  the  physical  process  associated  with  the  gaseous  source  had  a 
fixed  state  operator,  the  state  estimator,  through  its  finite  dimensional  implementation,  lend  itself  in  state- 
dependent  switching.  This  was  made  possible  through  different  evaluations  of  the  advection-diffusion  op¬ 
erator  over  different  grids.  Multiple  grids  representing  different  coarser/refined  grids  of  the  2D  spatial 
domain,  representing  different  regions  of  interest  within  the  2D  spatial  domain  were  used  to  evaluate  the 
advection-diffusion  operator.  The  resulting  switched  system  switched  to  a  different  state  matrix  depending 
on  the  current  location  of  the  sensor.  Thus,  a  performance  based  switching  that  bridged  computational  fluid 
dynamics  and  controls  was  considered  in  this  research  effort. 
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